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Abstract 

We present an exact analytical solution to a one-dimensional model of the Susceptible-Infected- 
Recovered (SIR) epidemic type, with infection rates dependent on nearest-neighbor occupations. 
We use a quantum mechanical approach, transforming the master equation via a quantum spin 
operator formulation. We calculate exactly the time-dependent density of infected, recovered and 
susceptible populations for random initial conditions, and compare our results with a low connec- 
tivity SIR model reported by Schiitz et al. 8)]. Our results compare well to those of previous work, 
validating the model as a useful tool for additional and extended studies in this important area. 
Our model also provides exact solutions for the n-point correlation functions, and can be extended 
to more complex epidemic type models. 
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INTRODUCTION 



The study of cooperative evolution of multi-agent systems is a part of many sciences, 
from biology and social science to physics, chemistry and engineering. The methods of 
statistical physics are often employed, and simple models are building blocks in this quest 
to understand complexity. From the point of view of non-equilibrium statistical physics, 
biology is an exciting area of investigation. After all, every living organism is an example of 
a far-from-equilibrium system, and stochastic processes are ubiquitous in biological systems. 

Epidemic-type models abound in the literature [lj], j^J, from very simple ones that capture 
the basic rules of the infection mechanism, to very complex models that account for spatial 
spread, age structure and the possibility of immunization It is interesting to see how 
some of these epidemic-type models have been applied successfully in other fields as well, 
such as social sciences (voter models, rumor spreading models) |4[,|5| or computer science 
(the spread of a virus in a computer network) 6|. Some models are deterministic, following 
a set of evolution equations with given initial conditions solved using the mean field theory 
approach, while others are stochastic and studied using methods such as the Langevin equa- 
tion, the Fokker-Planck equation and computer simulations. Some recent numerical studies 
of epidemic-type models can be found in |7|. 

Despite numerous studies and approximation schemes, exact solutions for epidemic 
models are rare. A study that sparked our interest was published in 2008 by Schiitz et 
al. jsj. This presents an exact solution for a stochastic one-dimensional SIR (suscepti- 
ble/infected/recovered) epidemic model. The method used is a quantum mechanical formu- 
lation of the master equation in terms of second quantized operators. The authors define 
cluster functions that describe the behavior of susceptibles adjacent to infected individuals 
at the cluster boundaries. They derive and solve exactily a set of coupled evolution equa- 
tions for these functions. Their exact solution shows the significant difference between low 
connectivity and high connectivity SIR models, and the role of fluctuations. Fluctuations 
are built into the exact solutions, but are missing in the mean field approach. 

Working towards an exact solution for an SIR model with higher connectivity, we study 
a variation of a one- dimensional SIR model in which we define different rates of infection 
depending on the number of infected neighbors. A susceptible individual with two infected 
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neighbors will have a different probability of being infected than a susceptible neighboring 
only one infected person. As in the traditional SIR model, we assume the possibility of 
recovery. (To differentiate it from other one-dimensional SIR models, we refer to our model 
as the dual neighbor model.) This model also can be cast as a non- conservative voter model, 
with the three classes of individuals defined as S - undecided, I - biased, R - decided. 

We here present an exact solution of the dual-neighbor SIR model. We employ a quantum 
mechanical approach to the problem, using the cluster function method used by Schiitz 
js]. The steady-state solution depends on initial populations of susceptible and infected 
individuals. It has fluctuations built-in, and has a stationary state different from the low 
connectivity SIR model. Although the overall trend of the solution is similar to that of the 
low connectivity model, there are significant differences as well. 

In Section 2, we define our model and its quantum mechanical representation. Next, 
we present the cluster function method and derive the evolution equations for the cluster 
functions and particle densities that fully resolve the model (Section 3). We conclude with 
an analysis of our solutions, summarize of our work, and suggest some interesting open 
questions (Section 4). 



DUAL NEIGHBOR MODEL AND ITS QUANTUM MECHANICAL REPRESEN- 
TATION 

The traditional SIR model consists of a fixed number of individuals N split into three 
classes: susceptible, infected, and recovered. In a high connectivity (mean field) model, each 
"node" of the network is represented by an individual in one of the three classes, in contact 
with every other node. A susceptible becomes infected with rate /3 when in contact with 
an infected; an infected individual recovers spontaneously with rate a; recovered individuals 
cannot change. At a particular time the average number of individuals in each class is 
represented by S, I, and R, with S + I + R = N. 
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The time evolution of these classes is governed by a set of coupled differential equations: 

dt 

dR - 

lu =+aI 

The first equation describes reduction of S via infection; the second shows I increasing via 
infection of susceptibles, and decreasing by spontaneous recovery; and the third shows the 
increase of R due to spontaneous recovery. This system of coupled nonlinear differential 
equations can be solved numerically, yielding the time dependence of each class of individ- 
uals, but without any information regarding correlations between individual nodes. This 
model exhibits smooth time dependence of the class populations, without any statistical 
fluctuations. 



In contrast, we propose a stochastic one-dimensional model in which each node is in 
contact with only two other nodes (as if arrayed along a line), and the rate of infection 
depends on the number of infected neighbors. Representing linear sequences of neighboring 
individuals via strings of symbols (e.g. SIS representing a susceptible to the left of an 
infected to the left of another susceptible), the dynamics of this model is defined as: 

ISI -»■ III with rate (3 
ISS — > IIS with rate A 
SSI — > SII with rate A 
I — > R with rate a 

The first process describes the mechanism of infection of a susceptible neighboring two 
infecteds; the next two processes describe infection with a different rate when the susceptible 
has only one infected neighbor; and the last process describes spontaneous recovery with 
yet another rate. We can pick the values of /3, A and a based on the application of the 
model. To model actual epidemics, for example, we would assign a higher rate of infection 
when a susceptible is in contact with two infecteds compared to just one, thus > A. 
Note that our model does not allow a succeptible individual with a recovered on one side 
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and an infected on the other (ISR or RSI) to become infected. As a model of disease 
transmission this represents an immunity provided by a recovered neighbor; in a voter model, 
it suggests that an undecided voter remains so if flanked by one biased voter and one who 
is decided. Practically, it says that the rate of infection in ISR and RSI configurations are 
negligible compared to other infection rates a and (3. Making such an assumption produces a 
Hamiltonian amenable to exact solution of the rate equations. As we shall exhibit, features 
of the results of this model validate it as a reasonable first approximation to real-world 
situations. 

The related mean field (deterministic) model is governed by the following differential 
equations for the average populations: 

^ = -pPs - MS 2 
dt_ 

^- = {3PS + \IS 2 -aI 
at 

dR - 

lu =+aI 

Compared with the traditional SIR model, the dynamics of this system is governed by 
three-point interactions between the S and I type individuals. This system also can be 
solved numerically for the time dependence of S, I and R. 



Proceeding beyond the mean field approximation, the time evolution of our model is best 
described by the master equation, expressing conservation of probability within a continuous- 
time dynamics. We let C represent a configuration, giving the state (S, I, R) of each of the 
N individuals. Proximity is defined by labeling each individual with an index % = 1, 2, . . . , N 
and assuming individual N is adjacent to individual 1 (periodic boundary conditions.) The 
master equation expresses the rate of change of the probability P(C,t) of finding the sys- 
tem in configuration C at time t as the rate of transfer of probability into C from other 
configurations less the rate at which C passes probability into others 9j: 

The transition rate r [C — > C'\ is the probability per unit time that configuration C changes 
into a different configuration C. 



Utilizing Dirac notation, we represent a configuration as |C), and use standard methods 



111 ] to build a vector representation for a general state as a probabilistic superposition of 



all configurations of a system: 

\P(t)) = J2 P (C,t)\C) (2) 

c 

where P(C, t) is the probability that the system will be found in configuration C at time t. 
This allows the master equation to be re-written as 

j t \P(t)) = -H\P{t)) (3) 

where the pseudo-Hamiltonian H has matrix elements: 

(C'\H\C) = -r{C -»• C), C'^C 

(C\H\C}= ^r(C^C'). (4) 

A formal solution to Eq. [3] can be written as \P(t)) = e~ Ht \P(0)) . In this formalism, the 
expectation value (at time t) of a physical quantity that has value M(C) for configuration 
C is 

< M > = p ( c > t)M(C) = < s\M\P(t) >, (5) 

C's 

where M = ^ c M(C)\C >< C\ is the operator corresponding to the observable M and the 
state \s > is the sum (with weight 1) of all configurations. Time dependence of < M > 
obeys 

9< Qt > =<a\[M,H\\P{t)> (6) 
where [M, H] is the commutator of the M operator with the Hamiltonian. 

The operators for our model can be written in terms of creation and annihilation oper- 
ators, following the rules for constructing a quantum Hamiltonian as presented in We 
define the a], to be the creation and annihilation operators for a type S particle at site 
i, and b\, bi to be the creation and annihilation operators for a type / particle at site i. A 
site with neither an infected nor a susceptible is assumed to contain a recovered. Operators 
at different sites commute, and operators of the same species at the same site anticommute. 
The number operator = a\di [Bi = b\bi) can only take values and 1, and its expectation 
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value equals the density of type A (B) particles. In this notation, 



- H =/3$>Ja t - A(l - Bi)}Bi_iB i+1 + A - A t (l - B^A^B^ + 

i i 

Xj^ibhi - M 1 - Bi)]Bi-iAi +1 + a ~ Bi]. (7) 

i i 

We can rescale the time variable in such a way that the rate A (assumed non-zero) becomes 
one. With that assumption, and defining 7 = f and 5 = f , the Hamiltonian obeys 

-H = J2[bhi ~ M 1 - B i)} [lBi-xB i+1 + A^B i+1 + Bi_iA i+1 ] + S^h ~ (8) 

i i 

EXACT SOLUTION: DERIVATION AND ANALYSIS 

We seek to find the time dependence of the average particle density of each kind of 
individual. To this end, we introduce n-point cluster functions following the method of |8|, 

K r (n) = < A r A r+1 . .A r+n - X B r+n >, (9) 
G r {n) = < B r - X A r ..A r+n - X B r+n > . (10) 

The Hamiltonian function (Eq. |8]) is invariant relative to translation along the string of 
sites, so if we assume a translationally invariant initial state, all future configurations also 
have this property. Translation invariance makes these clusters independent of their starting 
point r. Because a susceptible cannot change into an infected unless they are neighbors, 
clusters change only at their edges. 

The equations for the cluster functions are derived by calculating the respective commu- 
tators of the cluster operator products with the pseudo-Hamiltonian H. 
For n — 1, this gives 

= < s\[A r B r+1 ,H]\P(t) > = -SK(1) -tG(1), (11) 



dt 
dG{\) 

dt 



< s\[B r ^A r B r+1 ,H}\P(t) > = _( 7 + 2£)G(l) + 2G(2), (12) 



and for n > 1 

dK(n) 

dt " 
dG(n) 



dt 



< s\[A r . . . A r+n ^B r+n , H]\P(t) > = K(n + l)-(l + 5)K(n) - G(n), (13) 

< s\[B r ^A r . . . A r+n -xB r+n ,H]\P(t) > = 2G(n + 1) - 2(1 + 5)G(n). (14) 



The equations of motion for the particle density of susceptibles and infecteds are found 
in a similar fashion: 

d > = < s\[A r ,H]\P(t) > = -2K(2)- 7 G(1), (15) 
- < f r > = < s\[B r ,H]\P(t) > = 2K{2)+jG(l) -5 < B r > . (16) 

In order to solve for the particle density of the susceptibles and infecteds, we thus need first 
to find the cluster functions. 



Solution for cluster functions 

Introducing an operator s defined by the property sG n = G n+ i we can rewrite Eq{14] (for 
n > 1) as 

dG(n) 
dt 

which has a formal solution 



-2(1 + 5- s)G(n), 



G{n) = exp (-2(1 + 5 - s)t) G(n) t=0 . 

Representing the initial density of susceptibles and infecteds as rjs and r]j , we describe an 
initial state such that G(n) t =o = VsVl- From this, 

G{n) = exp (-2(1 + 5 - s)t) ^ = exp (-2(1 + 5 - Vs )t) V n s vl (17) 

Eq. [T2] for G(l)can be rewritten as 

+ ( 7 + 2 tf)G(l) = 2 exp (-2(5 + m )t) 4 v l 

and can be integrated by standard means to yield 

nsvj prised* -2 + 7 ) e-^t 
G{1) ~ 7^2^ • (18) 
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Using a similar method, we find the solution for K{n) to be 

K{n) = VsVie-r [1 - m r(l - e"*)], (19) 
where in order to express the solution more compactly, we have defined a relaxation time 



S+l-Tjs ' 



Particle densities 



The pieces are in place now to find the time dependence of the average number of suscep- 
tibles and infecteds. We assume that at t = there are no recovered individuals, therefore 
Vs + Vi — 1- Combining Eq.'s UHl [T8l and [T9l we find 

dt = [1 - mr{l - e~)\ - 7 -^-L-[2n 5 e^- + ( 7 - 2)e^] (20) 

where we have defined a second relaxation time r' = 1/(7 + 25). Integration yields 

^ /t ^ * 1 o 2 r 2 — - 1 2 2 / 7 \ — — 1 2717 — ^) 1 — * /oi\ 

< A > = ?7 5 + 2t] s r] I 5T e -+v S vM T + 7T~ ) e r + W// r e T ( 21 ) 

7 ~ ^/ 7 — *Vi 

where the integration constant 77^. turns out to be the stationary value of the density of 
susceptibles for large times. It is evaluated by enforcing < A > t =o= Vs- 

Vs = Vs l-Vi(vsr 2 (25 + r ]l )+r ]l — L_( T ^ + (7 - 2)^)1 . (22) 

This is non-zero because a susceptible can be "trapped" with a recovered on either side, and 
be no longer susceptible to infection. For the density of the infecteds, we combine Equations 
QH QJ and M to get: 

d<B> . „ 9 r t „ 9 9/ 7 \ a* 97(7 — 2) t . . 

dt + 6 < B > = Irf-^re-r + 2 V 2 sV j(r + ^-^)e~ + W ^_ 2 e' 7 . ( 2 3) 

leading to 



2r ?I^f /_ , 7 v -» „ 9 7(7-2) 



< £ > = C e - dt _ 2nS<Sre™ - — ^-ii-( r + ! )e"~ - 77577^7 — w ' r , e~^ (24) 

/5 5 + 277/ 7 - 2r// ; /5// (7-2r // )( 7 + 5) 1 ; 

where C is an integration constant to be evaluated using < B >t=o= 

C = VI + 245r + 2M-{t + — ) + W r tV + AV ( 25 ) 
5 + 2/7/ 7-2?7/ (7 - 2?7/)(7 + d) 

Note that the steady-state value for < B > is zero, due to the process of spontaneous 
transformation of an infected into a recovered, regardless of its neighbors. 

The time dependence of the average number of recovereds can be easily found from 

<i?> = l-<A>-<73>. 
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Analysis of solutions 



In an actual epidemic, major concerns are the time evolution of the number of infected 
individuals, and their maximum number. Basic questions include: at what time is the peak 
of the infection reached; what are the factors that control the maximum number of infecteds; 
and how long does it take for a population to fully recover? With these questions in mind, 
we examine the exact solution for the infected population, Eq. EU The solution depends 
on the initial fraction of infecteds, 777, the parameter 7 representing the relative rate of 
infection with two infected neighbors vs. that for one infected neighbor, and 5 representing 
the relative rate of spontaneous recovery vs. the rate of infection with a single infected 
neighbor. 

We point out two values of 7 that deserve special attention, 7 = 1 and 7 = 2. Physically, 
7 = 1 means that the the probability of infection is the same regardless of the number of 
infected neighbors. This most closely matches the previous SIR model js|. The case 7 = 2 
corresponds to an infection rate proportional to the number of infected neighbors, perhaps 
most closely aligned with the case of a more realistic three-dimensional medical infection 
model. 

Two typical time behaviors of the density of infecteds are shown in Fig. 1. Figure 1(a) 
shows the density rising from an initially small value to a maximum, and thereafter relaxing 
exponentially via spontaneous recovery until all infected individuals have disappeared. If 
the initial infection rate is high and/or the spontaneous recovery rate is high compared to 
infection rates, however, the model can exhibit time behavior that shows no peak, as in Fig. 
1(b). 

Figures 2(a) and 2(b) show how infection peak height and peak time vary with the initial 
state infection density. The peak height shows an increase as rjj increases, leveling off as it 
must as one approaches an initial state with the majority of the population already infected. 
The smaller the initially infected population, the longer delay until the maximum number 
of infecteds is experienced. For r/j values close to one, the peak will be early, if it appears 
at all. Because we have rescaled time in such a way that the infection rate of a susceptible 
with a single infected neighbor is 1, the proportional delay of the peak as a function of this 
single neighbor infection rate is hidden in these graphs. 
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Figure 3 shows how the spontaneous recovery rate 5 controls whether a peak in < B > 
appears or not. As expected, for a high spontaneous recovery rate the infected population 
begins to decay immediately once recovery is "turned on," as might be the case at the 
beginning of an immunization effort in the population. As the recovery rate drops, the 
figure shows the appearance of a peak that grows in size and in delay time. 

Figure 4 shows the comparison of our solution (continuous curve) with results of the 
low connectivity model obtained by Schutz et al. (dashed curve) and the mean field ap- 
proximation (dotted curve) for two representative sets of parameters. There are noticeable 
differences among the three solutions when 7 becomes large in comparison to 5 (Fig. 4(a)). 
These differences disappear for 7 values on the order of or smaller than 8, as shown in Fig. 
4(b). We notice that the peak of the infection in the mean field approximation is higher and 
happens at a later time than the one predicted by our model. 

The various behaviors exhibited in the figures well represent effects seen in actual situa- 
tions of disease spread [ij]] . This serves to commend the model as a reasonable approach to 
modeling such behavior, and the exact analytical solutions enable precise predictions about 
features of the solutions. 



CONCLUSIONS 

In this paper we analyzed a modified one dimensional stochastic SIR model, the dual 
neighbor model, with infection rates dependent on the nearest-neighbor occupation. Using 
a quantum mechanical approach and the cluster function method introduced in Schutz et 
al. js], we found an exact solution for the mean densities of susceptibles, infecteds and 
recovereds as a function of time and initial conditions. The quantum mechanical approach 
is a powerful analytical tool, because it leads to exact solutions not only for the mean number 
of infecteds, susceptibles and recovereds, but also for n-point correlation functions. 

We analyzed our solution for the density of infecteds in various parametric regimes, and 
also compared our solution with the low connectivity model reported in The qualitative 
behaviors exhibited by this model are straightforward and not unexpected. On the other 
hand, the details that this solution presents could lead to proposals for moderating the 
effects of an epidemic. This suggests the utility of further investigation of similar models, 
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perhaps in two dimensions. 

Although this model was introduced as an epidemic model, it can also be extended for 
other areas of study as well, such as voter problems, computer virus dynamics, or surface 
deposition. It can also be generalized to include correlated initial conditions. 

We hope to further our study to find exact solutions for more complicated (and realistic) 
SIR models that include time delay of infection, the possibility of immunization, and also the 
reappearance of the disease (R can become S). Unfortunately the cluster function method 
fails for this latter (SIRS) model, and other mathematical avenues must be pursued. 
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FIG. 1. Typical time evolutions of the density of infected individuals < B > for: (a) r/j = 0.1, 
Vs = 0.9, 5 = 1 and 7 = 2; (b) r?/ = 0.5, r? s = 0.5, 5 = 3 and 7 = 2. 
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FIG. 2. (a) Height of the infection peak as a function of the initial density of infected for 7 = 2, 
5 = 0.1; (b) Time of the infection peak as a function of the initial density of infected for 7 = 2, 
5 = 0.1. 
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FIG. 3. Density of infected individuals as a function of time and recovery rate 5 for r]j = 0.1, 
r]s = 0.9 and 7 = 2. 



16 




FIG. 4. Comparison with the mean field model (dotted curve) and the low connectivity model 
(dashed curve) jsj]: (a) For high infection rate 7 = 20 and r\j = 0.5, r]s = 0.5 and S = 2 ;(b) For 
low infection rate 7 = 3 and r\i = 0.5, tjs = 0.5 and 5 = 2. 
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